A knitting algorithm for calculating Green functions in quantum systems 
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We propose a fast and versatile algorithm to calculate local and transport properties such as 
conductance, shot noise, local density of state or local currents in mesoscopic quantum systems. 
Within the non equilibrium Green function formalism, we generalize the recursive Green function 
technique to tackle multiterminal devices with arbitrary geometries. We apply our method to 
analyze two recent experiments: an electronic Mach-Zehnder interferometer in a 2D gas and a Hall 
bar made of graphene nanoribbons in quantum Hall regime. In the latter case, we find that the 
Landau edge state pinned to the Dirac point gets diluted upon increasing carrier density. 



The field of quantum transport at the nanome- 
ter scale now includes a large number of systems in- 
volving very different physics. Examples include for 
instance mesoscopic devices in two dimensional het- 
erostructures [H, graphene nanoribbons [5], supercon- 
ducting weak links molecular electronic devices Q 
or ferromagnetic multilayers nanopillars Q . Although 
those systems have different structures and geometries, 
they are all quantum systems connected to the macro- 
scopic world through electrodes and, consequently, for- 
malisms developed to describe one of them can often be 
adapted to the others. This is in particular the case of 
the widely used Landauer-Biittiker formalism 0] which 
focuses on the scattering properties of the system. The 
formalism is very intuitive and general. However, it is not 
well suited when one is interested in what happens inside 
the sample or for performing a microscopic calculation 
for a given device. An alternative mathematically equiv- 
alent approach is referred as the NEGF (non equilibrium 
Green function) formalism [1, [^. NEGF, which is de- 
rived from the Keldysh formalism , provides a simple 
route to compute the physical observables from a micro- 
scopic model. It is now an extremely popular numerical 
approach to a very wide class of physical problems (see 
references in [IH). For instance, all the references men- 
tioned in examples above correspond to calculation done 
with this technique either from ab-initio or from phe- 
nomenological models P, i, i, II i, i, [H, [11 . At the core 
of NEGF is the calculation of the retarded Green func- 
tion G of the mesoscopic region in presence of the (semi- 
infinite) electrodes. A straightforward method consists of 
direct inversion of the Hamiltonian H. However, when 
doing so, one is restricted to rather small systems of a few 
thousand sites: for a system of size L in dimension d, the 
computing time scales as L^'^ while the needed memory 
scales as L^'^. An alternative algorithm, known as the re- 
cursive Green function technique [l|, d, , takes advan- 
tage of the structure of H to reduce drastically the com- 
puting time down to L^'^~^, putting systems of a few mil- 
lion sites within reach. In its original version (l^.[T5llT^. 
only the transport properties of the device could be com- 
puted, but recent progresses made it possible to get ac- 
cess to observables inside the sample (like local electronic 



density or local currents [H, [l3|) at a cost i^"^^^ in mem- 
ory. The recursive Green function technique suffers how- 
ever from a serious limitation: in its original formulation 
it is intrinsically one dimensional and most applications 
are done for quasi-one dimensional bars connected to two 
electrodes. On the other hand, real devices often have 
more than two electrodes and more complicated geome- 
tries. This paper is devoted to a general versatile algo- 
rithm to simulate multiterminal systems with arbitrary 
geometries and topologies. Earlier works in this direc- 
tion are scarce. Baranger and al. [ll considered the Hall 
effect in a 2D cross (a specific code was developed to 
handle this geometry). Modular algorithms [l^ [20] al- 
low to compute the properties of 2D quantum ballistic 
billiards. Most references with multiterminals involve di- 
rect inversions with small system sizes; others adaptation 
of the algorithm to the specific problem at hand [2l|, [24] . 
Competitors to the present algorithm are also being de- 
veloped [23, [2^ using alternative techniques such as dec- 
imation [25j . 

In this paper, we show that the recursive Green func- 
tion algorithm can be generalized to deal with arbitrary 
geometry, topology, number of connected electrodes and 
inner degrees of freedom (like spin for ferromagnets or 
electron/hole in superconductors). Our algorithm is con- 
ceptually simpler than the original as it is not based on 
a specific geometry. The sites are added one by one 
in a manner reminiscent of the knitting of a sweater. 
The algorithm is optimum in term of speed and a sig- 
nificant gain in memory is also achieved. The method 
allowed us to study the electronic Mach-Zehnder inter- 
ferometer [2^ [131 and anomalous quantum Hall effect in 
a graphene Hall bar \2§\ which have been the subject of 
recent experiments. Going from the first system to the 
second required virtually no additional development. 



I. NEGF FORMALISM IN A NUTSHELL. 

We consider a quantum system of N sites connected 
to several conducting electrodes. We use a general tight- 
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biding Hamiltonian for the system, 

N N 

where c| (q) denotes the usual creation (annihilation) 
operator of an electron on site i. The site index i stands 
for the position in space but can also include possible 
other degrees of freedom like spin or electron/hole. Uj 
is usually very sparse as only nearest (or next nearest) 
neighboring sites are connected. The electrodes / are 
semi-infinite systems at equilibrium with temperature T; 
and chemical potential /zj. They can be "integrated" 
out of the equations of motion and only appear in the 
formalism through self-energies that provide bound- 
ary conditions at the connected sites. We use an algo- 
rithm introduced by Ando [l^ [2^ for the calculation 
of those self energies. The physical observables can be 
simply related to the non equilibrium lesser Green func- 
tion Gfj{E) — i J dte~^^* {cjCi{t)) . For instance the local 
density pi and current reads, 

P« = ^Ini J dEG<{E), (2) 

/OO 
dE[U,G<{E)-tJ^Gf,m■ (3) 
-OO 

The simplicity of NEGF comes from the fact that, in 
absence of electronic correlation, G*- (which describes the 
system out of equilibrium) has a very simple expression 
in term of the retarded Green function G : = GS^G^ 
where I]< = Y.i M^l " ^i) and fi = 1/(1 + exp[(£; - 
lii)/kTi]) is the Fermi function. G itself is simply defined 
in term of the one-body Hamiltonian matrix Hij ~ tij + 
eiSij by, 

G{E) = l/iE-H-J2j:,). (4) 

II. KNITTING, SEWING AND UNKNITTING 
ALGORITHMS 

A. Basic idea behind recursive algorithms. 

The problem of computing G{E) is thus reduced to 
a conceptually simple task, finding the inverse of the H 
matrix (to which we implicitly include the self energies). 
Our basic tool is to judiciously divide H = Hq + V be- 
tween an unperturbed part Hq and a perturbation V. 
Typically, the perturbation will be the hopping elements 
that allow to glue two separated part of the system to- 
gether. The Dyson equation G — g + gVG which relates 
G to the known g — 1/{E — Hq) is the corner stone of all 
recursive algorithms. Suppose that (i) one is interested 
only in Gap for a small subset of sites labelled by a, /3 
(the electrode sites for instance) and (ii) Vij only con- 
nects a small number of sites labelled by The follow- 
ing glueing sequence allows one to get the needed matrix 



elements in three steps. (I) Dyson equation restricted to 
the connected sites (via Vij) is a close equation: 

Gy = gij + X gikVkiGij (5) 

kl 

The number of those connected sites being small, the G^ 
can hence be easily computed. (II) In a second step one 
obtains the elements 

Gaj = gaj + X gakVklGlj (6) 
kl 

and (III) in the third step 

Gq/3 = gaf} + ^'^gakVklGip (7) 
kl 

are computed. In the original recursive Green function 
algorithm the above sequence is used in the follow- 
ing way: one considers a bar of width M and length 
L. The bar is sliced in L stacks and the perturbation 
Vij are the hopping elements that connect the different 
stacks. The system is then built recursively as the stacks 
are added one at a time and at each step the glueing 
sequence is used. The gij are known either from the pre- 
vious calculation (for the bar side) or from a numerical 
direct calculation for the added stack. 



B. Knitting algorithm for global transport 
properties. 

Our knitting algorithm is based on the same glueing 
sequence, but the sites are added one by one. We start 
by indexing the sites according to the order in which they 
are going to be added to the system. Fig. [T] shows a car- 
toon of a typical system together with the notations used 
in the following. The main difficulty of the algorithm lies 
in the book keeping of the various Green function ele- 
ments and precise definitions are compulsory. At a given 
stage of the knitting, when we have already included sites 
1 ... ^ — 1 and are about to include site A, we distinguish 
between four categories of sites labelled by different in- 
dices: 

• the connected sites are the sites that are directly 
connected to A via hopping elements i^o- • They are 
labelled by the index a and appear in very small 
number (less then the number of neighbors of a 
given site A). 

• The interface sites labelled by i,j are the sites that 
still miss some of their neighbors, i.e. that will 
be themselves connected sites later in the knitting. 
Hence, the Green function elements for those sites 
is to be kept in memory. Note that it is a dynam- 
ical definition, i.e. at each step, some sites appear 
and others disappear from the interface. The to- 
tal number Af of interface sites scales as a surface 
M cx L''-!. 
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FIG. 1: The system scheme and notations. The sites are la- 
belled according to the order in which they are added to the 
system. The boxes define the three electrodes coupled to the 
mesoscopic system. The various letters stand for the added 
site (A), electrode sites (a,/3, 7), connected sites (ai,a2) and 
interface sites as discussed in the text. Bold circles indi- 
cate sites whose GF elements are updated at current knitting 
step. The thick dashed line separates the part already in- 
cluded (left) from the part that is still to be knitted to the 
system (right). 

• Updated sites, the sites whose Green function ele- 
ments are updated at each step of knitting. They 
belong either to the interface or to the electrodes. 
They are noted by a, (3; bold circles on Fig. [TJ 

With these notations, we can apply the glueing sequence 
and express the Green function with the added site 
in term of the Green function G'-^^^I of the system com- 
posed of ^ — 1 sites. (I) The first step reads, 

G[^i = 1/iE -^A-J2 tA.d^-%A) (8) 

aa 

Note that it is the only place where we actually perform 
an inversion and it is done on a scalar quantity. (II) The 
second step reads, 

41 = T.y^^-''t^^GAl (9) 

^[4/3 ^ ^^AA'^AaG^fg ^\ (10) 

(III) The last step concentrates almost all the computing 
time (oc M^), 

^[^1 _ ^[A-i] ^[A] 1 ^[A] 

^AA 

Note that the previous formula has a very simple physi- 
cal interpretation in term of paths: the amplitude for an 
electron to go from site /3 to site a is the amplitude which 
avoids site A plus the amplitude which goes through site 
A. The factor 1/g|^Jj removes the double counting of 
the loops from A to itself. Once the glueing sequence 
is completed, the interface is updated: the new site is 
added while sites that now have all their neighbors can 
be removed. For instance in Fig. [U once site A has been 
added, site ai can be deleted from the interface. The 



previous sequence is iterated until all the N sites have 
been added to the system. Eventually we get G^^ of the 
entire system. However, we only get the matrix elements 
between the (few) sites a, (3 connected to the electrodes. 
Those matrix elements give access to all transport prop- 
erties like conductance or shot noise, but no informa- 
tion on what happens inside the sample. The computing 
time (needed memory) of the knitting algorithm scales 
as M'^N (M^) in agreement with the original recursive 
algorithm for a wire. 



C. A sewing algorithm for calculating local 
observables. 

We now proceed with extending the knitting algorithm 
to the calculations of local observables. In practice, we 
need matrix elements of the type Gar between electrodes 
and any inner site of the system r. Such an extension 
has been done in [l| for the original recursive algorithm: 
one performs a first recursive calculation and saves the 
partial Green functions Gar for future use. When the 
calculation is completed, one starts a new recursive cal- 
culation, beginning from the other hand of the system. 
Along the way, one recovers the saved Green functions 
for the left part of the system and uses the glueing se- 
quence to glue them with the (freshly calculated) Green 
function of the right part of the system. This scheme can 
be generalized to an arbitrary geometry as well: we first 
perform a full knitting calculation. Then we start back- 
ward the sewing algorithm and sew the sites A one by 
one in reversed order (i.e. from N to 1). We label sites 
with an index smaller than A (" left" part of the system) 
by indices without prime {%,],...) and the sites (strictly) 
higher than ^ ("right") by prime indices ■■■)■ Intro- 

ducing the "interface self energy" Stj = J2i'j' tij'Gj'i'ti'j, 
the glueing sequence reads, 

GaA = + ^[^^ ^^jG^^ , (12) 

Ga'A = Ga'j'tj'jClJ . (13) 

j'i 

The result of the calculation is then used to update Sij 
and one can proceed with A— 1 and so on. The drawback 
of this algorithm is that many (~ NM) matrix elements 

[rl 

Gar must be stored during the first knitting calculation 
which limit practical calculations to a few hundred thou- 
sand sites. 



D. Saving memory with unknitting. 

The last piece of algorithm, called unknitting, allows 
to recalculate the matrix elements G^ in the backward 
calculation instead of saving them. Indeed, using Dyson 
equation for Hq = H — V, it is possible to "remove" sites 
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from the system and express G^-^ ^1 as a function of G^-^l 
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AA 



(14) 



The above equation should however be taken with care. 
Indeed the interface is not the same for Gl-^-i] and GI-^1 
so that some of the matrix elements G^^\ on the right 
hand side are not stored in memory. In the bulk of the 
system the interface is of constant size so that there is 
only one site Ti- that belongs to the interface of G^-^"^' 
but not of G^'^l. The matrix elements for this site can 
also be recomputed, 

G^an ^' = G^;^\{tTZAG AaY^ - X! '^^aa^^^'t'yAtnA 

which completes the algorithm. Theoretically, the un- 
knitting algorithm allows to decrease the memory from 
MN to max{N, A'P), hence completely removing the 
memory bottleneck in the calculations. We found that 
it is indeed the case in the middle of the tight-biding 
band where all the channels are conducting. Outside this 
region, the last equation above is however numerically 
unstable and introduces an error in the calculation that 
increases as exp(aiV/M). The origin of this instability 
can be understood from the example of a simple perfect 
ID chain with unit hopping. In this case, Eq.® takes 
a simple form GI-^I = - GI-^-^I). When \E\ > 2, 

i.e. the chain is no longer propagative (only evanescent 
waves can be found) and this equation converges toward 
a simple attractive fixed point. After several iterations, 
the convergence is achieved up to numerical precision and 
the initial condition is lost: the equation can no longer 
be inverted to follow our steps back. We found that the 
unknitting algorithm is nevertheless useful. Depending 
on the precision needed, the matrix elements g[^ can 
be saved in the forward calculation only every few sites 
instead of every site. In practice, in the worth case (bot- 
tom of the band where almost all modes are evanescent), 
we found that a factor 10 in memory could be gained 
while keeping the numerical precision better than 10"^". 



III. APPLICATION: MACH-ZEHNDER 
INTERFEROMETER IN A TWO DIMENSIONAL 
ELECTRON GAS. 

As a first application of our algorithm, we consider a 
Mach-Zehnder interferometer similar to the ones stud- 
ied in recent experiments ,27] . The devices are made 
in a two-dimensional high mobility GaAs/AlGaAs het- 
erostructures which we characterize by its mobility /i, 
electronic density Us , size C and perpendicular magnetic 
field B. The simulations are performed by discretizing 
the Schrodinger equation on a square grid with a lattice 
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FIG. 2: Mach-Zehnder interferometer in a 2D gas modelled 
by a scalar tight-biding model. Mobility is = 5 - lO^cm^/Vs, 
electron density — lO^^cm"^ and magnetic field B = 0.2 
T. (a); local current intensity when a bias voltage is applied to 
lead and the other contacts are grounded (1.2 million sites, 
blue colors corresponds to no current and red to maximum 
current), (b): differential conductance dIz/dVo as a function 
of the number of flux quanta <j!> through the hole, (c): visibility 
coefficients Vi and V2 (see text) as a function of the electron 
mobility /i. 



spacing b. The resulting tight-biding model has near- 
est neighbor hoppings t and disordered on-site poten- 
tial which are chosen randomly within [— 1^/2, +W^/2]. 
The magnetic field is added within a discretized Lan- 
dau gauge by adding a phase $ in the hoppings ele- 
ments X along the y direction: t^, = ie*^'^*"=" 

(where {nx,ny) is the position on the grid along x and 
y). The energy E is measured from the bottom of the 
band, and we found no significant deviations from the 
continuum limit as long as we kept E < O.lt. The 
tight-biding parameters are related to the experimental 
ones as follows: = E/{2ntb'^), B = <^h/{eb'^) and 
^ = (e//i)967r&^(t/W^)^. The latter formula was obtained 
from a calculation of the Drude conductance of our tight- 
biding model and was checked by direct numerical calcu- 
lations. The total number M — L^ris of electrons in the 
sample and the conductance per square g = efins are 
both independent of the discretized step b. The middle 
of the n*'' plateau of quantum Hall effect is found for 
n = E/{4TTt<S>). 

The Mach-Zehnder interferometer is the electronic 
analog of the well-known optical device. The system con- 
sists of a loop connected to four contacts, one of which lies 
in the center part of the loop, see Fig. [5^. The physics 
involved is fairly straightforward: the system is placed 
under a high magnetic field in the quantum Hall regime 
at the first Hall plateau (the current is supported by one 
edge state). All contacts are grounded except contact 
which is placed at a slightly higher voltage Vq. The in- 
jected current follows the edge channel until it reaches 
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FIG. 3: Same as Fig[2^ for two dirtier samples with mobility 
/i = 3 • lO'^cm^/Vs (a) and ^ = 2.5 ■ lO^cm^/Vs (b), electron 
density Ub = lO^^cm"^ and magnetic field B = 0.2 T. Zoom of 
(a) shows the current that escapes the first contact, and makes 
the second loop. It is responsible for the presence of a second 
visibility harmonic Vi = 1 V2 = 0.1. Similar zoom of (b) 
shows backscattered current which is recovered by contact 2. 
Vi = 0.99 V2 = 0.01. 



a first QPC (Quantum Point Contact which works as 
a perfect beam splitter) and is split into two parts, see 
Fig. [5^. The two edge states are eventually recombined 
at the second QPC and the current J3 is collected in 
contact 3. Along the way, the two edge states pick up 
a difference of phases that includes the magnetic flux 
through the hole. Note that Fig. [^h- is not a cartoon, 
but an actual calculation of the local density of current 
injected from lead 0. Along the way, the two edge chan- 
nels pick up a phase difference which produces interfer- 
ences. Fig.[5]D shows the differential conductance dl^/dVo 
in unit of e^//i as a function of flux ip (in unit of h/e) 
along with a (almost perfect) sinusoidal fit. Our data 
are in complete agreement with what is obtained from 
Landauer-Biittiker theory. In particular, the visibility Vi 
decreases when the transmission T of one QPC departs 
from 1/2 as Vi = 2^T(1 - T). 

To proceed further, we increase the disorder in our 
sample and measure the visibility of the interference 
pattern as a function of the mobility fj,. We find that 
dl^ / dVo{(j)) is well fitted by including two harmonics, 

dis /dVo cx 1 + Vi cos 2tt4) + V2 cos 47r0 (15) 

The parameters Vi and V2 are shown in Fig. [2t for a 
typical sample. For > 4 • lO^cm^V^^s"^, there is no 
backscattering in the sample and the interference pattern 
is perfect Vi = 1 and V2 = 0. Below 4 • lO^cm^V^^s^^ 
however, backscattering sets in and the visibility de- 
creases strongly. More importantly, the visibility be- 
comes extremely sensitive to the disorder configuration 
and huge sample to sample fluctuations are observed. Si- 
multaneously, the second harmonic V2 sets up. Those 
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FIG. 4; Quantum Hall effect in graphene. Hall conductance 
a^y and longitudinal resistance p^x are plotted as a function 
of inverse magnetic field 1/$ (a) and carrier density n (b) in 
presence of a small disordered potential (10% of the hopping 
matrix elements). In (a) carrier density is 0.18 electrons per 
atom while in (b) a magnetic flux "1> = 0.014/i/e per hexagon 
is applied. The inset of (b) shows a zoom of the transition 
between plateaus, (c) shows the local current intensity when 
current is injected from both contact 1 and 4 (3- 10* hexagons, 
iV = 2) 



results are readily understood by looking at the current 
intensities injected from lead 0, as shown in Fig. [3] for 
two different samples with mobilities slightly lower then 
in Fig. [2^. When the disorder is strong enough, the edge 
channels get significantly disturbed and are able to reach 
(partially) the opposite wall (as seen in Fig. [5]d). This 
backscattering leads to a decrease of the visibility. At the 
same time, a path that avoids contact 1 appears which 
causes V2 to increase (the corresponding path makes 
a complete tour of the central island, see the inset of 
Fig. [3^). It is notable that no second harmonic V2 was 
observed experimentally which is, as the samples had mo- 
bihty of the order of 10^cm^V~^s~^ or higher [2^, \2T\ . 
consistent with our findings. We conclude that static 
disorder (or more trivially, a bad central ohmic contact) 
cannot be at the origin of the reduced visibility in these 
experiments. 



IV. APPLICATION: ANOMALOUS QUANTUM 
HALL EFFECT IN A GRAPHENE HALL BAR. 

In this second application, we consider a Hall bar in 
graphene in the quantum Hall regime, see Fig. |4l:. We 
used the standard tight binding approach ^0] for the 
graphene hexagonal lattice, where each site represents a 
carbon atom and is connected to his neighbors by a hop- 
ping t. We use zigzag edges for the system and the leads 
which are semi infinite nanoribbons. To avoid reflection 
at the lead/system interface, the leads also include the 
magnetic field. The magnetic field is included within 
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FIG. 5: Comparison of the edge channels for a graphene 
(a) and a semi-conductor heterostructure (b) two-dimensional 
electron gases. The plots show the local current density in the 
y direction dIy(x)/dV as a function of the distance x from the 
border of the sample for four consecutive Hall plateaus corre- 
sponding to A'^ = 0, 1, 2, 3 for graphene (a) and n = 1, 2, 3, 4 
for the heterostructure (b). All graphs are translated for clar- 
ity from bottom to top. Magnetic field B = 14 T. 



V. CONCLUSION 

We have presented a set of algorithms that allows to 
calculate the transport and local properties of a generic 
class of tight-biding models. Using the unknitting tech- 
nique, some significant gain of memory could be gained 
with respect to previous techniques. However, our main 
point here is the global simplicity of the algorithm. We 
could easily implement it to get a versatile code that 
treats complicated geometries, such as the ones presented 
in the applications, with the same ease as one treats the 
usual quasi-one dimensional bar, ubiquitous in the liter- 
ature. 



APPENDIX A: IMPLEMENTATION TIPS. 



a gauge similar to the one used for the Mach-Zehnder 
above, by introducing a phase $ in the hopping elements 
between the A and B sites. Noting a (1.42A) the dis- 
tance between carbon atoms, the electronic density Us 
and magnetic field B are related to E/t and the total 
flux $ (in unit of h/e) per hexagon as follows: = 
{E/atf^l{9TT) and B = 2h^ / {?,^/?,ea^) . The middle of 
TV*'* Hall plateau is found for N + 1/2 = E"^ / {2^/i^l^) 

Fig. 3] shows the transverse conductance and longitu- 
dinal resistance as a function of carrier density (a) or 
inverse magnetic flux (b) for a graphene Hall bar. Our re- 
sult for a mesoscopic Hall bar are numerical counterparts 
to the experimental data of [28| . In particular we recover 
the presence of plateaus in the transverse conductance at 
quantized values {N + 1/2) • Ae^/h. The = plateau 
corresponds to a Landau level pinned to the Dirac point, 
and is special to the Dirac equation symmetry class. We 
also observe mesoscopic fluctuations at the transition be- 
tween plateaus similar to those known in usual 2D gases, 
see inset of Fig. 0^. The edge current density at iV = 2 
injected simultaneously from contact 1 and 4 is plotted 
in Fig. Uj: which is the graphene counterpart to Fig. [^t. 

We note one important difference between the edge 
states of graphene and those of a regular heterostructure: 
while Fig.HJ; has been made for N = 2, i.e. for the third 
plateau in the sequence, one observes only two peaks in 
the local current density, i.e. apparently two edge states 
only. This is revealed more quantitatively in FigIS] where 
we have plotted the cross section of the current density 
for the first plateaus of both the graphene {N = 0, 1, 2, 3) 
and the heterostructure (n = 1,2,3,4). The area under 
those curves correspond respectively to 1/2, 3/2, 5/2 and 
7/2 in the graphene and 1, 2, 3 and 4 in the heterostruc- 
ture. The number of peaks however is respectively 1, 1, 
2 and 3 and 1, 2, 3 and 4 so that the iV = peak (corre- 
sponding to the E = Landau level) gets blurred behind 
the other peaks upon increasing carrier density. This 
peculiar behavior could possibly be observed in STM or 
local compressibility experiments (3l| . 



In this appendix, we have grouped a few technical 
points which can help to get significant gains in com- 
puting time, or maybe more importantly in human de- 
velopment time. 

1. The computing time for adding one site in the knit- 
ting algorithm is dominated by the last step of the glueing 
sequence, Eg. pTjl and scales as O(M^). Hence, almost 
all the computing time of the algorithm is concentrated 
in a single line of code which can be optimized aggres- 
sively. The best results were obtained by using low level 
BLAS routines for this simple step. Note that the only 
inversions that take place in the algorithm are scalars: at 
no point we need to rely on matrix inversion routines in 
contrast to the original recursive technique. 

2. A practical difficulty is to maintain dynamically the 
list of active interface sites. We use the following simple 
algorithm: to each interface site is associated a list of 
the site's neighbors that have not yet been knitted to the 
system. As new sites are knitted to the system, this list 
is updated. When a site has all his neighbors knitted, he 
can be removed from the interface. 

3. The algorithms described in this paper allow to 
solve a very generic class of tight-biding models. Almost 
as important as those solvers is the ability to easily con- 
struct new systems. In our implementation, we represent 
internally a system as a generic graph: each site possesses 
the list of its neighbors, together with the corresponding 
matrix elements. A set of routines is used to operate 
on those graphes. A very important one is a function 
that allows to stick two systems together by identifying 
some sites of the two systems. This routine allows us 
to construct a system from small pieces in a "Lego" like 
way. For instance, the system of Fig. [^t is obtained in 
four lines of code by "sticking" four rectangular systems 
together. 

4. The intensive pieces of the code should be written 
in a low level language (in our case CH — h). However, 
we have found that very significant gain in developing 
time could be obtained by making use of the code in 
a higher level scripting language. In our case, we use 
Python, and the automatic wrapping of the CH — h code 
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into the Python API has been done with SWIG. The de France" as well as from the EC Contracts IST-033749 
systems constructions, calculations and plotting are done "DynaMax" is acknowledged, 
in small Python programs that constitute our input files. 
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